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This work presents a novel approach to describe spectral properties of graphene layers with well 
defined edges. We microscopically analyze the boundary problem for the continuous Bogoliubov-de 
Gennes-Dirac (BdGD) equations and derive the Green functions for normal and superconducting 
graphene layers. Importing the idea used in tight-binding (TB) models of a microscopic hopping that 
couples different regions, we are able to set up and solve an algebraic Dyson's equation describing 
a graphene-superconductor junction. For this coupled system we analytically derive the Green 
functions and use them to calculate the local density of states and the spatial variation of the 
induced pairing correlations in the normal region. Signatures of specular Andreev reflections are 
identified. 

PACS numbers: 73.23.-b, 74.45. -l-c, 74.78. Na, 73.20.-r 



o 

> 

ON 



X 



I. INTRODUCTION 

Ideally, an interface between a graphene layer and a su- 
perconductor provides a unique opportunity to study the 
interplay between a relativistic-like band-structure and 
pairing interactions between electrons and holes. The 
advances in the techniques of isolation and manipulation 
of graphene layers together with the possibility to deposit 
nanoscale superconducting electrodes on top of them is 
bringing this idealized situation closer [l| . 

Several theoretical studies of the electronic and trans- 
port properties of graphene-superconductor junctions 
have been reported in recent years. As in the case of 
normal metals, a key ingredient for understanding these 
properties is the concept of Andreev reflection [2[ by 
means of which an incident electron on the interface is 
reflected as a hole while a Cooper pair is created on the 
superconductor. In graphene, however, Andreev reflec- 
tions may couple electrons in the conduction band with 
holes in the valence band, which leads to a specular re- 
flection, as opposed to the usual retro Andreev reflection 
0. Very recently we have shown that interface bound 
states (IBSs) appear for energies within the supercon- 
ducting gap due to the interplay of Andreev and normal 
reflection in graphene-superconductor junctions 

Within these studies various different methods and 
approximations have been applied. In particular, the 
Bogoliubov-de Gennes-Dirac equations (BdGD) have 
been solved under the approximation of neglecting in- 
tervalley scattering and applied to the calculation of 
the differential conductance and spectral properties in 
graphene-superconductor (GS) junctions [3|, |5|. Other 
works have extended this method to situations includ- 
ing an insulator (GIS junctions @) or to the study of 
the Josephson effect in SGS junctions 0- In spite of 
their interest these works are limited in the sense that 
they cannot provide an answer on the effect of the atomic 
scale properties of the interfaces or what would happen 



in the presence of well defined edges in the graphene lay- 
ers. This is an important issue because it is well known 
that the electronic spectrum of a graphene nanoribbon 
is strongly dependent on the direction in which its edges 
are defined. Armchair and zig-zag cases have been the 
most studied but one can also envisage and study edges 
along arbitrary directions. From the point of view of 
experiments it is also becoming possible to study trans- 
port properties in graphene nanoribbons with well de- 
fined edges at the atomic scale Q. 

These type of studies typically require the use of mi- 
croscopic tight-binding (TB) models and some of them 
have been conducted in the case of GS junctions 1,1131. 
Although such methods are usually implemented numer- 
ically, much progress has been made in the derivation of 
analytical expressions within TB models combined with 
Green function techniques 0. However, these expres- 
sions are typically more complex than those which would 
be obtained within the BdGD model. 

The aim of the present work is to go one step beyond 
the existing works [llj and develop a method based on 
the BdGD model which could deal with the presence 
of well defined edges or interfaces at the atomic scale. 
Our method starts from the determination of the Green 
functions for the BdGD equations in the case of iso- 
lated graphene or superconducting layers. In order to 
fix the boundary conditions for these isolated regions it 
is in general necessary to mix solutions corresponding 
to different valleys. For the purpose of illustrating the 
method we would consider separately the cases of arm- 
chair and zig-zag edges. For obtaining the properties of a 
coupled graphene-graphenc or graphene-superconductor 
region we define a short range coupling potential which 
mimics the usual hopping term in the TB models and 
then solve the corresponding Dyson equation analytically 
for several junctions. As a test we show how known re- 
sults are recovered in different limits. We further use 
the method to study the spatial properties of the local 



density of states (LDOS) and the induced pairing cor- 
relations on a graphene layer. We show that, on top 
of an atomic scale modulation, these spatial correlations 
exhibit signatures of specular Andreev reflections. They 
also display a long range decay due to the presence of the 
above mentioned IBS. 

This work is organized as follows; in the next section 
we implement Green's functions for graphene nanorib- 
bons with armchair or zigzag edges. In section III it is 
explained how we model the microscopic link between 
layers and solve the corresponding Dyson's equation. Fi- 
nally, section IV covers the coupling between a normal 
sheet of graphene and a superconducting one. Analyt- 
ical expressions for the Green functions of the coupled 
system are given. Profiles for the LDOS and the spatial 
behavior of the induced superconducting pairing correla- 
tions are presented as well. This paper is closed by some 
concluding remarks. Furthermore, an appendix has been 
added to show details of the calculations. 



II. GREEN'S FUNCTION APPROACH FOR 
GRAPHENE LAYERS WITH EDGES 
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FIG. 1: Geometry of the systems treated in this work: a hor- 
izontally finite sheet of graphene with armchair edges (left) 
or zigzag edges (right). The sheets are assumed to be infinite 
in the vertical direction. The edges are located at positions 
XL and Xr. These positions can be taken to be infinite, rep- 
resenting a bulk of graphene (infinite plane) or a semi-infinite 
layer (infinite half-plane). For each geometry we represent the 
chosen Dirac points in the Brillouin zone. 



Graphene is a single layer of carbon atoms arranged 
in a honeycomb structure. Its hexagonal lattice is con- 
structed as the superposition of two triangular lattices 
A and B. A nearest neighbor TB model for graphene 
can be, in the continuum limit, represented by a lin- 
earized Hamiltonian with six Dirac points at each cor- 
ner of the hexagonal Brillouin zone. Only two of these 
points are inequivalcnt which, for the orientation of the 
layer shown in the left panel of Fig. [TJ can be assumed to 
lie in the x-axis and denoted as K± = {zLK,0). Single- 
particle excitations around these points are represented 
by a Dirac Hamiltonian like H± = hva± ■ k — Ep&o = 
hv [±kax + QO'y] — EpCTQ, where dx^y arc the Pauli ma- 
trices (with ctq the 2x2 identity matrix) acting on 
graphene's lattice subspacc, Ep is the Fermi energy and 
the velocity v = y/^tga/2h « 10^ m/s is proportional to 
the lattice constant a = 0.246 nm and to the nearest- 
neighbor hopping energy « 3 eV on the honeycomb 
lattice of carbon atoms. With this choice K = Air/Sa. 
The degree of freedom due to the A~ B lattices define a 
pseudospin which points in the direction of motion of the 
particle and the pseudospin operator a± = (±0-3;, CTj,) is 
used to define the Hamiltonian on each valley. Finally, k 
is the wave vector measured from each Dirac point K±. 

Wave functions around each Dirac point satisfy the 
equation H±^± {x,y) = E^± {x,y), with E > being 
the excitation energy of an electron-like quasiparticle. In 
the present work we consider clean samples which are 
infinite along the y-direction but may have edges or in- 
terfaces on x-direction (see Fig. [1]). Thus, momentum 
along the j/-axis (hq) is conserved and the wave function 
can be written as ^±{x,y) — e'^'^y (j>±{x) . With the re- 
placement k(g) = {—idx,q) in Dirac's equation, linear 



independent solutions for each valley are 
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and e^'" = hv(k ± iq)l{Ep -\- E). The constants cf 2 are 
determined by boundary conditions for the wave func- 
tions (/)± {x). 

A complete description of the low energy electron ex- 
citations on a graphene layer is reached by superposing 
solutions on both valleys. Each solution $± is modulated 
by a rapidly varying plane wave from each valley e'^± '' 
[1^ thus giving 

$ (x, y) = e'««0(x) = e^^y [e'^"0+ (x) + e-*^"0_ (x)] . (4) 

This wave function is solution of Schrddinger's equation 
for a Hamiltonian that combines both valleys. However, 
it is helpful to use a four-component spinor notation for 
the wave function made from the two-dimensional spinors 
for each valley in which the phase factors are omitted (i.e. 
7/)^ = {(1)1, 0^) = (0^, 0f , 0^*, 0^?)) [il. This new nota- 
tion is useful in order to calculate the full Green function 
of the system, including both valley and pseudospin de- 
grees of freedom. The four-dimensional spinor obeys the 
equation Hijj (x) = Eip (x), where 



is the Dirac Hamiltonian in sublattice and valley spaces. 
We are using the notation • " • for 4 x 4 matrices and • " • 
for 2 X 2 matrices. 

The Green function associated to = e^'^yip{x) 

is (x, x' , y) = J (x, x'; q) e^'^^dq, where G^ {x, x'; q) 
satisfies the 4x4 matrix equation 



H [q) - EI] G^ {x, x'; q) = 5{x - x')i 








-•S{x-x')i, 



(6) 



with / the 4-dimensional identity matrix. Hereafter we 
implicitly assume that E stands for E + ir] and that the 
limit 77 — > has been taken to obtain the retarded com- 
ponent of the Green function. From the elements of the 
full Green function G^ we can define a valley superposed 
Green function G^, associated with the wave function of 
Eq. (HI), as 



(7) 



G^{x,x\y,y')^ j dqe'«(^-^') x 
E e<'^^-'^-''-'^')G;''{x,x';q). 
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Besides allowing valley mixing, this Green function in- 
clude the phase factors expi (K^ • r — Kj/ • r') which are 
crucial to describe the presence of well defined edges or 
interfaces at the atomic scale. This valley superposed 
Green function is derived in the next sections for arm- 
chair and zigzag edges. 

The method that we implement to calculate the Green 
function is based in using the asymptotic solutions of the 
differential equation ([6]) . Basically, we extend the method 
used in Ref. to the relativistic Bogoliubov-de Gennes 
equations. Thus, the full Green function in pscudospin 
and valley spaces can be written as 



G^ (x, x'j q) = 
E _A^,ij'^{x)-r/{x')--f ■,x<x' 

E' A'^.^'^ix)-r/ix')-j -xyx' 



(8) 



where the indexes /i and v represent different valley 
projections of the asymptotic solutions. We build the 
Green's function using tensor products as G{x ^ x') oc 

^<(>)(''') ' '/'>"('<) (^')- Where ■0<(a;) and i/'>(2;) are the 
asymptotic solutions that fulfill the boundary conditions 
to the left and right edges of the system, respectively. 
In Fig. [2] (a) the possible processes for an armchair 
edge are shown. The labels and v are always deter- 
mined by the valley index of the incident particle. On the 
other hand, ip'^f^yj{x') represent asymptotic solutions for 
the transposed Dirac's equation with the same bound- 
ary conditions and opposite wave number ■(/'^^^^(r) ~ 

exp(— jg2/)V'<(>)(a;). Since H'^ {—q) = H{q), the wave 

functions '!/'<^>)(a;) fulfill the same differential equation 
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FIG. 2: (Color online) Schematic representation of the asymp- 
totic solutions of Eq. ([6]) for a system with armchair edges. 
Scattering at an armchair edge always changes the valley pro- 
jection. The processes with non zero probability combine so- 
lutions on the left edge with one valley projection (e.g. -0^) 
with solutions on the right with the opposite valley index 
(i/;>). (a) Example of valley mixing at armchair edges {xl 
and Xr) for incident particles with a fixed valley index, where 
solid (dashed) lines represent particles with valley index -|- 
(— ). Reflection processes for electron (b) and hole-like (c) 
quasiparticles. For a semi-infinite region, xr goes to infinity 
and the asymptotic solutions are outgoing plane waves. For 
zigzag edges we have the same behavior without the change 
in the valley index. 



than i/jjij^j (x). The constants A^'^ are determined from 
the condition 



E A'^,^^{x)-r/{x)-j 



(9) 



hv 



which is obtained integrating Eq. ([6|) around x = x' . 

In our definition of the Green functions in Eq. ([8]) we 
have introduced the parity matrix 7 in order to make the 
scalar product ^Gip invariant under Lorentz transforma- 
tions and spatial inversion (parity transformation). Here 
■0 is the adjoint Dirac spinor defined as V' = The 
transformation rule for Dirac spinors is -0 — > "0' = Sip^ 
where 5 is a 4 x 4 matrix that fulfills S^^S = 7. With 
this rule, the product (^ipip) = 'ipS'^^Sip' = 'tjjip is Lorentz 
invariant (for a complete description of the symmetries of 
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Dirac fermions see Ref. [T^). Now, since G oc ^ip'^j, we 
have that G ^ G' = SG-^S-^ and the product (i/'Gi/') ' = 
Sj (^SG^Sj^ jtp — ip'^^Gip = 'ipGip is also Lorentz in- 
variant. The matrix 7 in reciprocal space changes both 
valley and pseudospin indexes. Then, since the spinor is 
in the basis B^, the matrix 7 can be writ- 

ten as 7 = Tx^cFx and we can check that the Hamiltonian 
H is transformed with 7 as "jH (p) j = H (— p). In the 
next section, where the armchair edge is studied, we use 
this form of the matrix 7. A simplified form of matrix 7 
for one valley is derived in Appendix [X] and is used for 
zigzag edges. 



A. Armchair edges 

We can now apply this method to derive the Green 
functions for a layer of graphene with armchair edges. 
Within the geometry depicted in the left panel of Fig. 
[U the direction perpendicular to the edge may be infi- 
nite or have a finite size W. In this later case, we have 
two edges at positions xl and xr. The main character- 
istic of an armchair edge is that scattering at the edges 
changes the valley index. This is a direct result of van- 
ishing either the wave function (Eq. Q) or its derivative 
on both lattices at the armchair edge. Asymptotic solu- 
tions for this system correspond to incoming waves from 
one valley that are refiected on the other valley. All pos- 
sible processes are illustrated in Fig. [D^a). The index 
± symbolizes an incoming wave around the Dirac point 
K-t that is refiected at the interface. For incident waves 
refiected at the left (<) or at the right (>) edge we have 
the four-dimensional spinors in valley and lattice spaces 

V'<(»(a;) = 

( T ) + ^ii-f'-'' ( ) (10) 

±ikx ( Vl(2) 





^1(2) ) ^ '^(«)'^ V 



(11) 



The refiection amplitudes '"^(^j are calculated impos- 
ing boundary conditions to the two-dimensional wave 
function of Eq. (g]) at the edge x^^f^), this is = 
4>{xr) = or dx4>{xL) = dx(l){xft) = 0, which yields 



'r ~ 



(12) 



where s = ± corresponds to the case in which the wave 
function (— ) or its derivative is vanished at the in- 
terface. The main differences between both cases are 
studied in the next section. For simplicity we have cho- 
sen that Xl S (—00,0] and x-^ G [0,oo). Combining 
the asymptotic solutions, as has been explained in the 
previous section, we obtain a general expression for the 



full Green function with armchair edges (G^™) equiva- 
lent to Eq. ([8]). The constants A'iXi' in this equation 
correspond to processes in which the excitation changes 
its valley projection when it propagates from one edge 
to the other. Therefore, these processes have zero prob- 
ability and, from Eq. one explicitly obtains = 
= 0, A+_ = A'_^_ = -i/ {2hvcosa{l-r+r]^)) 
and = A'_^ = -i/ (2hv cos a (l - rlr^)). Sub- 

stituting these results in G'^"^ we can derive the valley 
superposed Green function for a layer of graphene with 
armchair edges 



2hv cos a (1 — Tj^) 



Gr'{x,x') 

_)_j,-^+g-«(^f+fc)h-2;'| _^ ^-^i(K+k)(x+x') 



i{K+k)\x-x'\ 



_j^r+^~i(K+k)(x+x') 



2hv cos a (1 — ^iffi) 



' ^-t(K-k)\x-x'\ _|_ j.+ j.^f,KK'-k)\x~x'\ 
^^+^'^iK-k)ix+x') ^ ^-^^^(K-k)(x+x')^ ^^^t ^ (^g) 

The poles of the Green function are determined by 
cosa = 0, (l-r£r+) = and (l - r+r^) = 0. The 
first one results in the bulk dispersion relation hv\q\ = 
\Ef + E\. In the second and third cases we reach the 
condition 2{K± + kn)W — 2mT, with W = xr ~ xl and 
n an integer. Thus, the allowed values of the transverse 
momentum for a finite armchair sheet are /c„ = ^ T |f ; 
independently of the boundary condition chosen for the 
armchair edges. This result is in agreement with previ- 
ous works on graphene nanoribbons with armchair edges 

We now consider the case of a superconducting 
graphene region, where the description of electron and 
hole excitations is done within Nambu space. We assume 
a s-wave pairing which leads to a constant gap A which 
is diagonal in sublattice space. The system is then repre- 
sented by the BdGD equation written in lattice, Nambu 
and valley spaces (8x8 matrix). It reads 



A* 




















= E 




\4>\) 











H--E% 
A* 





A 

Ef-H^ 



(14) 



with E positive unless otherwise specified. Each valley 
Hamiltonian is written as H± = —ihv[-^dxCrx-\-dyay\. 

Two-dimensional spinors (jf^ represent electron or hole 
like excitations for each valley in lattice space. Whenever 
the pairing potential is assumed constant, the low energy 

spectrum is given hy E — J A'^ + (j^p — hv^Jk^ + . 
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We define the transversal momentum as hvk^ ,^ 



{hvq) , with n = V£^2 _ The so- 
lutions of the BdGD equation are a direct product of 
graphene's bulk solutions and the usual BCS solutions. 
The pair potential couples electron and hole-like excita- 
tions with opposite momentum, which means different 
valley index, thus (j)^ and (j)^ correspond to different val- 
leys. This allows us to reduce the degrees of freedom of 
the problem to Nambu and pscudospin spaces. 

In order to study a semi-infinite superconducting 
graphene region with an armchair edge extending from 
XL = to Xfl — >■ oo, we define a finite and constant pair 
potential A(a;) = 6(a::)A. Reflection amplitudes are then 



fixed to 



r and 



0. The 



asymptotic solutions are 



u[v) 
v(u) 



ie{h) 



I Q 



u{v) 
v(u) 





V2e(l/i; 



(15) 
(16) 





Vle(2h) 



u{v) 
v{u) 



(x) 



(17) 
(18) 



where ipie{h), with 
([3|) making a — > 



V2e(l/i)/ \V{U) 

= 1,2, is obtained from Eq. 
'^e{h) ^^'^ defining e'""''') = 



hv (^fcf(^-) + iq^ /\E^ ± n\ and (v^) = (1 ± n/E) /2. 

The function ip'^'^'^ describes an incoming electron-like 
(hole- like) quasiparticle with energy E > Ep [E < Ep) 
and valley index ± that is reflected at xl into an electron- 
like (hole-like) quasiparticle with valley index =p. 

The 8x8 matrix Green's function can be written as 
the superposition of all possible quasi-electron and quasi- 
hole injection processes depicted in Figs. WO^) and^c). 
The resulting Green function, written in Nambu, valley 
and pseudospin spaces, has the form 



Gl, [x < x') = 

j=e,h 



< X' 



E A;^,^ir(x).Vi:'^^(x')(7®fo) ;x>x' ' 

j=e,h 

(19) 

where Tx,y,z and ro are Pauli matrices in Nambu space. 
The parity matrix 7 is diagonal in this space and only 
acts on graphene's valley and pseudospin degrees of free- 
dom. 

As in the normal case, the coefficients A'^i^^I with ^ = v 
vanish as they correspond to processes which change the 



valley polarization without reflection at the edges. Ad- 
ditionally, the processes that couple quasi-electron with 
quasi-hole and vice versa have zero probability. This 
can be explicitly found from the condition given by 
Eq. © where A^/}i^ = 0, ^;_(_+) = ^-+(+-) = 
iEe+(-^"^- 1 {2hvQ.cosae) and = ^+-(-+) = 

ii?e+(^''"''/ (2fi,fil cosa/i). Thus, after applying bound- 
ary conditions and solving Eq. ([9]), the valley superposed 
Green function in Nambu space for a superconducting re- 
gion with an armchair edge is 



Al-Al 



with the following definitions 



(20) 



Ae(h) = fe(h) {cos ^Q;f(^)CT0 +tanaf(^)(7j^) ±i5e(/i)'5-a; 

= cos[X(x'-x)]e±'^'(Ml-'-H 
+r cos [K [x' + x)]e±''==<'» (^'+^) 
<7e(h) - sgn(x'-x)sin[A"(x'-x)]e=^*'==('')h'-H 
+r sin [K (x' + x)]e^'<^) (^'+^) . 

It is interesting to take this expression to the heav- 
ily doped limit in which Ep ^ E^ A, hvq and thus 
cosaf^ — ?> 1 and sinaf^ 0. Within this approxima- 
tion, for a microscopic distance xq ~ a from the interface 
we have 

GI'""^^ {xo,xo) « ^ {[1 + rcos2A:xo] (To ^21) 
(n^o + ^Tx) + irsin2A'xoO-^ ® f^] ■ 

For the boundary condition of vanishing the derivative 
of the wave function (r = 1), we can evaluate the Green 
function at the edge (xq = 0). Thus, 



G^'™(0,0) 



2i 



{Eao (E)To + Ao-o (E) t^) , (22) 



which corresponds to the momentum independent BCS 
Green function without sublattice structure Q. Then, 
2/hv corresponds to an averaged Fermi energy DOS per 
unit length for the superconductor. 

On the other hand, for the boundary condition over 
the wave function (r = —1) the Green function vanishes 
at xo = 0. For finite xq, the Green function can ex- 
hibit a structure in pseudospin space given by the last 
term in Eq. (j2ip . This expression is used in the follow- 
ing sections when we microscopically couple normal and 
superconducting regions. 



B. Zig-zag edges 

The geometry of the graphene sheet is now set to have 
a zigzag edge along the y-axis. With this new orientation 



6 



of the layer the Brillouin zone rotates as illustrated in the 
upper part of the right panel of Fig. [Tl which allows to 
select the Dirac points at K± ~ {0,ztK). We still call 
hq the conserved momentum along the y direction. We 
consider either a finite layer with edges at positions xl 
and xr or a semi- infinite one in which xl or goes 
to infinity. Zigzag edges are formed by a line of atoms 
all pertaining to one sublattice. These edges do not mix 
valleys so we can treat them separately and use a Dirac 
Hamiltonian of just one valley, iJ+ = hv {kax + l^^y) ~ 
EpCTQ. Thus, the asymptotic wave functions including a 
reflection at the edge are 



(23) 
(24) 



The reflection amplitudes depend on the atoms chosen 
to form the zigzag edge. Thus, if we choose the border 
at Xl to be formed by atoms of A [B) lattice, the one at 
xji is formed by atoms pertaining to B {A) lattice (see 
right panel of Fig. [1]). Imposing the boundary conditions 
4'<{xl)\a(b) ~ '?^>(^fl)l_B(A) ~ we obtain the reflection 
amplitudes: 



' L{R) 



(25) 



When combining the asymptotic solutions to build the 
Green functions, as it was done in the previous section, 
the elements of the full Green function that mix valleys 
(H — , — h) are zero. In what follows the full Green func- 
tion is understood to represent only a projection over 
one valley (i.e. G^^'^^). Following the steps given in the 
previous sections we can write 



Gl^ix,x') 



A4i^{x) ■ (j)'^^{x') ■ ;x<x' 
A'0>(x) • (f>'^{x') -7 ■,x> x' 



When restricted to a single valley, the parity matrix 7 is 
equivalent to (t^ (see Appendix [X] for a detailed discus- 
sion). For a ribbon of graphcne with the left zigzag edge 
of type B (and thus the right one of type A) we have 



2hv cos ( 



ptk{x'-x)^^^-\^ ^Bj,Ap-tk{x'-x)^^^-\^_^_ 



X < 



^ik{x+x 



"flfl 



:x < X 



:x > X 



(27) 

For the other valley we obtain the same result with the 
change ipi,2 ^ ¥'2.1 • As it was explained in the previ- 
ous section, cos a = gives the bulk dispersion relation. 
Furthermore, if the layer has a finite width which is set 



to W 
into 



xr — Xl, the condition 1 



2ikW 



ik 
ik 



transforms 



(28) 



For k real, this expression leads to the quantization 
of transverse momentum in the ribbon (i.e. q = 
kn/ ta.nknW). On the other hand, if the transverse mo- 
mentum is a pure imaginary number (k ~ —iz) Eq. 
transforms into 

-2zW _ 1~ ^ 



V, (29) 

q + z 

whose solutions correspond to surface states, at E = Ep, 
localized along the edges of the ribbon [l^ • 

Analogously to the armchair case, we now consider a 
superconducting region spread over the x > infinite 
half-plane with a zigzag edge at x = 0. If we couple 
electronic excitations with K-|_ index to hole-like quasi- 
particles with the other valley index, Eq. (|14p is reduced 
to a 4 X 4 matrix equation in Nambu and pseudospin 
spaces 



A* 



E^ 



(30) 



The asymptotic wave functions for quasi-electron and 
quasi-hole injection are written in Nambu and pseudospin 
space as 



e(/i) 



(X) = { 



Ah) ±ikx,^ ^ ^ 
e 'fle{2h) 



}«(" 



u (v) 
(u) 



Alikx 

e Vle(2h) 



u (v) 
V (u) 



(31) 
(32) 



where the edge has been chosen to be formed by atoms 
from sublattice B. Boundary conditions for the wave 
function at the zigzag edge determine the reflection am- 
plitudes. Thus, for electronic excitations we substitute 



(26) a — af in Eq. ((25|) . while we change a 



-af for hole 



excitations. The resulting Green's function for a semi- 
infinite superconducting region reads 



where 



(33) 



^le{2h}^2e{lh) 



2 cos a'~ 



V2e{lh)V2e(lh) 
'Ple{2h)V>\e{2h) 



;x < X 
■,x>x' 



When Eq. ((33)) is evaluated at the edge of the graphene 
layer it reduces to 



Gf^(0,0) 



2hv 




(34) 



7 



Furthermore, in the heavily doped hmit {Ep ^ 
E, A, hvq), the previous expression reduces to 



Gf^(0,0) 



(35) 



where the sign identifies the cases with x < x' ox x > x' . 
Thus, in this Hmit we roughly have a BCS Green function 
projected to the site A. 



III. DYSON EQUATION FOR COUPLING 
SEPARATE REGIONS 



Once the expressions for bulk, armchair and zigzag 
graphene sheets, as well as for the semi-infinite super- 
conducting regions have been obtained, we concentrate in 
the determination of the Green functions for the coupled 
system. This amounts in principle to solve the integral 
Dyson equation 



G^{x,x') = g^{x,x') 

+ / dxidx2g4,{x,xi)V{xi,X2)G^{x2,x'), 



(36) 



where now denotes the valley superposed Green func- 
tion of the coupled system, corresponds to the uncou- 
pled ones and V is an appropriate perturbation describ- 
ing the coupling between the two regions. In general, we 
must use valley superposed Green functions in Dyson's 
equation to have a microscopic description of the cou- 
pling between graphene layers. In order to simplify the 
model while still keeping the possibility to describe inter- 
faces of arbitrary transparency we shall assume that the 
perturbation only connect points on an atomic scale close 
to the interface. We thus define a general perturbation 

V (Xi, X2) = PtgUeffS {Xi - .To) 5 {X2 + To) (Tj;, (37) 

where /3 G [0, 1] is a dimensionless parameter that con- 
trols the transparency of the interface and xq is a micro- 
scopic distance from the interface. The effective hopping 
matrix tga^ffax has a structure dictated by the change 
of sublattice associated with the hopping between the 
two graphene sheets at the interface. In this expression 
a.eff represents an "effective distance" which we deter- 
mine from the condition that the bulk graphene result 
is recovered when coupling two semi-infinite regions with 
/3=1. 

We first discuss the case of two semi-infinite graphene 
sheets with armchair edges. The simplest possible choice 
would be xo — > 0**", which can be applied when the con- 
dition of vanishing the derivative of the wave function at 
the edge has been taken. In this way one may convert 
Eq. p6|l into an algebraic equation, which for the local 
Green functions at the interface can be written as 



G 



R(L) 



where R, L denotes (x, x') 

by 



(0^, 0"*=) and gL{R.) are given 



9l{r) 



-2i 

hv 



cos ^ a&o + tanaiTy) 



(39) 



One can check that the bulk Green function is recov- 
ered for tteff = •\/3i/4, which therefore sets the condition 
of perfect transparency for this case. The reader is ad- 
dressed to Appendix IB] for more details about the bulk 
results for a normal or superconducting graphene region. 

A slightly more cumbersome matching condition has to 
be introduced if one would like to reproduce the limiting 
results of the TB model for such an interface Q. This 
requires using the Green functions obtained by vanish- 
ing the wave function at the edges of the graphene sheet 
and take xq = a/4. One obtains then a Dyson equation 
like dSll) with R, L denoting (a;, x') ^ (±f , ±|) and gL,R 
given by 



9l(r) 



-\/3i r 



2hv 



^/3 (cos ^ a&Q 



tan aa. 



V) 



(40) 



These expressions exactly reproduce the TB results in 
the continuum limit, this is a — >■ but with Ka — > 47r/3 
[1]. One can further check that for this case the bulk 
graphene result is recovered for Ce// = a/2 which is the 
distance between lines of atoms along the x direction in 
an armchair edge. 

We now focus on zigzag edges. In this case one can 
disregard the valley mixing but care should be taken due 
to the discontinuity of the Green functions aX x = x' . As 
discussed in the previous section we select the condition 
of vanishing one of the wave function components at the 
interface and take — 0"*". The precise value of ~ a 
is here irrelevant because the Dirac points lie along the 
y direction. 

When we evaluate Eq. ([27]) at the edge of each region 
we get 



9l 



9R 




z gj^ ;x < x' ^ 
E 5^ ;x' < X ^ 0^ 

= g^ ; X < x' 0^ 
= 5^ ; x' < a; — > 0^ 



(41) 



(42) 



(38) 



where we have implicitly assumed a A[B) termination 
for the L{R) graphene sheet. The Green function of the 
coupled system must be constructed including the discon- 
tinuity of the uncoupled Green's functions at the edges, 
i.e for x,x' < 

Gl {x, x') = g (x, x') + {tgUefff g {x, 0~) ct^, 
X 1^1 - (tgfle//) gjifTxglax^ 9r^x9{^ 

where ^(a;,0-) = cxp (— ifcx) g]!^, gL{^~,x') = 
exp (—ikx') (jj^ and g (a;, x') corresponds to Eq. (|27p with 
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r? = 0. When x < x' , and using that for this case 

'L 



-r^, this resuhs in 



Gl ix,x') 



2hv cos a 



(44) 



where we have replaced tg = 2hv/{-/ia). The bulk re- 
sult is recovered at perfect transparency (/? = !) for 
fle// — a/So/S, which corresponds to the horizontal sepa- 
ration between atoms of the same sublattice. This result 
corresponds to the projection over the point K_|_. For 
the other valley, the same result with tfi o ip2 is reached 
and the bulk result is recovered for perfect transparency 
superposing both valleys. 



IV. GRAPHENE-SUPERCONDUCTOR 
JUNCTION 

All the ingredients necessary to analyze the graphene- 
superconductor coupling have been already introduced. 
In the present section, we analytically derive the Green 
function of the coupled system both for armchair and 
zigzag edges. We locate the interface at a; = 0, parallel 
to the y axis, with the normal region covering the x < 
infinite half-plane. Since the superconducting region 
is represented in Nambu space, in order to use Dyson's 
equation as has been explained in section III, the Green's 
functions of the normal region arc written as 



9L 



9e 
9h 



(45) 



Electron excitations are states above the Fermi en- 
ergy. Then, corresponds to Eq. p3p for arm- 
chair edges or Eq. ((27|) for zigzag edges if we set 

Qfe = a, with hvke = \J {Ep -\- E)^ — (hvq)^ . On the 
other hand, hole excitations have energies below Ep. 
Thus, gh corresponds to the same equations but with 
the change ah = arcsm hvq / {Ep — E), with hvkh = 

'\J [Ep — E)"^ — [hvqf' . For the superconducting region 
(i?), Eq. ([20|) for armchair edges or Eq. ([33]) for zigzag 
edges, evaluated at the interface, stand for in Dyson's 
equation. 

The coupling with the superconductor appears in the 
Green function as a perturbation of the semi-infinite re- 
sult (jL given by Eqs. (fT5|) and ([?7|) . Thus we can write 



Gl{x, x') = x') + 5G{x, x'), 



(46) 



where 5G is the correction with respect to the uncou- 
pled case. Our formalism allows in principle to obtain 
analytical results for 5G in the case of arbitrary effective 
hopping and coordinates x and x' . 

For the case of an armchair interface, using the bound- 
ary conditions of vanishing the derivative of the wave 



function and setting /3 = 1, for x = x' we obtain 



5Gf-^{x^x') 
cos 2Kx 
-i sin {2Kx — ag) 
I Warm ( sinae 



hv COS Oe 

-i sin {2Kx + ae) 
cos 2Kx 



(47) 



—I 
sin ttp 



where we have defined the auxiliary quantities 



E 



A2 



2 sin 1 + TT cos a/i + —r (sin a^, — sin ah) 



02 



E 



Daim = 2 + 2cOSaeCOSQ;;i -f 2— (cOSQ!e -I- COSa/i) 

+ (1 + COS [ae + ah]) ■ 

The denominator Z^arm contains information on the 
spectral properties of the coupled system. In particu- 
lar, the condition Z^arm = is satisfied by states with 
energy given by 

E (e(^^+^'^)/2 - sign(£;2 - )e-(^^+^'-)/2) 

— ^ ±^ / ^ '-, (48) 

A 2vcoshAe cosh A/i 

where we define Ae.^ = sign(q)arcosh(?iu(7/|£' ± Ep\) for 
evanescent electron or hole states with \hvq\ > \Ep ± E\. 
As demonstrated in Ref. these states correspond to 
interface bound states (IBSs) arising from the interplay 
between normal and Andreev reflection at the junction 
region. The decay of these states into the normal region 
is given by e^/?=''', with £,e,h = hv/ (|£; ± £^^1 sinh Ae,h), 
which can be clearly much larger than the BCS super- 
conducting coherence length ^ = hv/A when Ep <^ A. 
In the opposite limit, Ep ^ A, the states are confined 
to the interface. 

On the other hand, for the calculation of the LDOS 
it is necessary to obtain the diagonal components of the 
local Green's function ([ee, AA] and [ee, BB]). Whenever 
Kx = TiTT, they simplify to 



'~'L\ee,AA 



{x,x) 



hv 



1^ 

I cos Oie 



'~'L\ee,BB 



{x,x) = 
tan Qfp I 



(49) 



Furthermore, for analyzing the induced pairing corre- 
lations in the normal graphene layer, the electron-hole 
elements of 6G arc needed. Within the same conditions 
of the previous result, but for any values of x and x' , 
these elements reduce to 



Gaim 
L\eh,AA{BB 



{x,x') 



hvD^ 
E 



{cos [K {x — x') 



1 + cos [og — ah] + — (cos ae + cos ah) 

1 . , COSQfe - COS ah 

+ — sm [K [x ~ X )\ 

2 cos ae cos ah 



E . . 

sm ae + sm ah + sm [ae 



ah 



(50) 
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We now analyze a normal-superconductor graphene 
junction with a zigzag edge at a; = 0. As it was explained 
at the beginning of this section, we have chosen the nor- 
mal (superconducting) region to be at x < (x > 0) 
and to have the edge formed by A (B) atoms. For these 
conditions we obtain the following correction 



hvD^, 



■0' 



2 e" 



COS dh 

O \ COS Cte COS Oih I ' ' COS O.^ COS Qh 



(51) 



with D,, = l + 



f3' 



Again, when this denominator is set to zero one obtains 
the dispersion relation for the IBSs along this edge. 



E 
A 



=.(Ae + A,0/2 



±- 



/32sign(£:2 _ Ej,)e 



-(Ae + A^)/2 



(52) 



which is in agreement with the result of Ref . [3| ■ 

On the other hand, the electron-hole component of the 
Green function is given by 



(53) 



A. LDOS for armchair and zigzag edges 



of Fig. [3] has been calculated for values |a;| /a = 3n, with 
n a positive integer. Similar profiles are reached for dif- 
ferent multiplicities of the coordinate x. The oscillatory 
behavior at the atomic scale is shown in the inset of this 
panel for the boundary conditions of vanishing the wave 
function. 

On the other hand, the left panel of Fig. [3] shows the 
LDOS for a zigzag edge at a distance inside the normal 
region d = —0.4^. The results correspond to different 
values of the parameter j3 controlling the interface trans- 
parency. For the uncoupled case with = 0, a localized 
edge state appears at i5 = Ep (dashed line in Fig. [4]). 
However, for a non-zero transparency the edge state splits 
and evolves into a band inside the superconducting gap. 
As it was demonstrated in Ref. [1], this band is formed 
by the IBSs which give rise to the peaks in the DOS for 
energies around E ~ when /3 = 1. 

It is important to emphasize that although the surface 
state a,t E = Ep for zigzag edges has been thoroughly 
studied, our results demonstrate that the superconduct- 
ing proximity effect splits this state and produces a shift 
that depends of the transparency between the normal 
and superconducting regions. 

The LDOS for the armchair case with the same set 
of parameters is illustrated in the right panel of Fig. H) 
As can be observed there is no localized state for the 
uncoupled case (/3 = 0) and when /? increases, two peaks 
arise around E = ±A^. 



As a test of the method we analyze in this section 
the LDOS for a graphene-superconductor junction. The 
electronic LDOS per unit length for a distance x = x' = d 
is defined as 



-1 

p{E,d) = — / dqTi-[lm{GL\,eiq,E;d,d)}] 

^ J —OO 



(54) 



The left panel of Fig. [3] shows the LDOS for an arm- 
chair edge close to the interface {d = —0.1^) obtained us- 
ing the boundary conditions of vanishing the wave func- 
tion or its derivative. For comparison wc also show in 
this figure the results obtained for the TB model calcula- 
tions of Ref. 1]. The resuhs of Fig. [3] are normalized to 
the density of a bulk graphene layer with zero doping at 
E = A, po = ^AjW-v^ (see Appendix |B]). The discrep- 
ancies between the different matching conditions tend to 
disappear when we move away from the interface, as it 
is shown in the right panel of Fig. [3l This panel also il- 
lustrates the oscillatory behavior of the LDOS inside the 
normal region of graphene around the bulk value (indi- 
cated by the dashed line). The period of the oscillation is 
given roughly by fmjE and the amplitude decreases with 
the distance to the interface. When studying the spatial 
evolution of the LDOS of an armchair layer of graphene, 
it is well known that there is an oscillatory behavior at 
the atomic scale [§, [l^. This is a direct result of the 
valley mixing that happens at an armchair edge. For the 
sake of clarity, the evolution of the LDOS with the direc- 
tion transversal to the interface shown in the right panel 



B. Induced pairing correlations 

The coupling with the superconductor induces cor- 
relations between electron and hole-like excitations in 
the normal region by the proximity effect. We can 
map the spatial variation of these correlations using the 
Green function of the coupled system. The element 
GL\e.h,AA(x^ x' , y) corrcsponds to the injection of electron 
excitations at the point x' on the y = axis and its 
propagation as a hole excitation to the rest of the plane 
(x, y) after an electron-hole conversion at the interface. 
The probability Pe^h that an electron injected at (x', 0) 
would be converted into a hole at {x, y) is proportional 
to 



d9e^'"G'L|e/.,ss' {q,E-x,x') 



(55) 



with S, S' = A, B. 

In our study of the spatial variation of the supercon- 
ducting correlations we set x' = — ^. We find it interest- 
ing to explore the change in the induced pair correlations 
for electrons injected within the gap [E < A) when the 
doping level varies from Ep ^ A (left panel of Fig. O to 
Ep ^ A (right panel). As it was demonstrated in Ref. 
|3|, these two limiting cases correspond to the regimes 
in which Andreev reflection exhibits respectively a retro 
or a specular character. In these limits and in the case 
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200 




FIG. 3: (Color online) LDOS in the normal region with Ep = for a semi-infinite layer of graphene with armchair edges 
coupled to a superconducting region (left panel) . Profiles for both boundary conditions (full blue line when vanishing the wave 
function and full red line for the derivative) are equivalent to the TB results (dashed line). The right panel shows the oscillating 
behavior of the LDOS for E = 2A along the direction transversal to the interface. The transversal distance d is normalized to 
the BCS superconducting coherence length ^. All the results show the same oscillating behavior around the bulk value. The 
oscillatory behavior at the atomic scale is shown in the inset for the boundary condition of vanishing the wave function. 




FIG. 4: (Color online) LDOS profiles, at a distance d = —0.41^ inside the normal region with Ep = 0, for different transparencies 
of the zigzag (left panel) and armchair (right panel) interface. In the zigzag case, the edge state splits and evolves into a band 
inside the superconducting gap when the transparency fi goes from (uncoupled case represented by a dashed black line) to 
1 (perfect transparency, solid black line). For the armchair case the superconducting proximity effect is manifested by the 
appearance of sharp peaks in the LDOS for energies _E ^ A. 



of perfect transparency, the main difference between an 
armchair edge and a zigzag one is that the atomic oscil- 
lations happen along the x ot y axis, respectively. Thus, 
we only analyze the results for an armchair edge with a 
fixed multiplicity. In Fig. [5] these two cases are shown 
for an incident energy E = 0.75A. In the case Ep > A, 
shown in the left panel, electron-hole conversion occurs 
mainly for ~ decaying fast for y ^ This confine- 



ment of the pairing correlations on the y direction is a 
consequence of the underlying diffraction pattern for the 
injected electrons at one point which has a spatial exten- 
sion of the order of hv/Ep. The behavior is drastically 
different in the limit Ep < E < A for which specular 
reflection is enhanced (right panel). In addition to the 
electron-hole conversion for normal incidence, the results 
exhibit the presence of long-range superconducting cor- 
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1 2 3 4 5 




1 2 3 4 5 




12 3 

FIG. 5: (Color online) Spatial mapping of induced pairing 
correlations in the normal region of graphene. The results 
show the squared absolute value of the Fourier transform of 
GL\eh,AA for au incident energy E = 0.75A and x' = The 
panel on the left shows the case Ep > A with Ep = 2.1A, in 
which retro-reflection is enhanced. The right panel shows the 
opposite case Ep < A, with Ep = O.IA, in which specular 
reflection dominates. 



relations along the y axis, provided that x ^ These 
long-range correlations can be directly associated with 
the presence of IBSs which in this regime have a spatial 
extension inside the normal region which is much larger 
than ^. At the same time, the extended correlations can 
be interpreted as a signature of divergent Andrecv reflec- 
tion trajectories characteristic of the regime Ep ^ A. 



V. CONCLUSIONS 

In this paper we have developed a method that allows 
to obtain the Green function in inhomogeneous graphene 
systems with well defined edges, such as nanoribbons and 
graphene-superconductor interfaces. In this approach the 
asymptotic solutions of the Bogoliubov-de Gennes-Dirac 
equations satisfying the appropriate boundary conditions 
are used to analytically build the Green functions asso- 
ciated to the quasiparticles in the system. 

Within this formalism we have analyzed the LDOS of 
a normal-superconducting graphene junction and showed 
how it is affected by the type of edge and the interface 
transparency. In particular, for the zigzag case, when 
the parameter /3 (which controls the transparency of the 
junction) increases, the localized edge state appearing at 
E = Ep splits and evolves into symmetric bands within 
the superconducting gap. These bands correspond to the 
interface bound states already described in a previous 
work (3]. 

Furthermore, we have studied the spatial distribution 
of the induced pairing correlations inside the normal re- 
gion. This correlations illustrate the crossover between 
retro and specular Andreev reflection. We indicate the 
appearance of long-range correlations between distant 
points on the graphene layer close to the interface due 
to the presence of IBSs. These superconducting correla- 
tions are enhanced in the regime corresponding to spec- 
ular Andreev reflection. 

We would like to conclude stressing that the analytical 
results obtained in the present work could be very use- 
ful for the study of transport properties in more complex 
hybrid systems which combine finite normal and super- 
conducting graphene regions. Work along these lines is 
under progress. 
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Appendix A: Green's functions for zigzag edges. 
One valley description. 

For the case of a graphene sheet with zigzag edges, 
we have set the Dirac points to be at K± = (0, zb/C). 
Thus, the Hamiltonian on each valley is given by = 
hv[kax^q(jy\ — Ep&o. The Hamiltonian of the full 
system is related with the one given in Eq. ([5]) by 
H^^ = THf, with 



The matrix 7 = f^; (g) (j^; is then transformed as %2 = 
T^T = Ty ® &y. When we apply this matrix to we 
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have 7zz (-&zz (k)) 7"^ = i?zz(— k), which corresponds to 
a parity transformation. 

Since the zigzag boundary conditions do not mix val- 
leys we can work only with one valley (K+) using Hamil- 
tonian H™ = + . For this Hamiltonian the effect of the 
parity transformation, restricted to the sublattice sub- 
space, is equivalent to set 7 = (Tz. 



Furthermore, the bulk LDOS is obtained integrating this 
simple result, thus giving phuik{,E) = 2{E + Ep) lf?v^ . 

On the other hand, vanishing the reflection coefficient 
(r) in Eq. (|20p we reach the bulk solution for an infi- 
nite superconducting region. The local valley superposed 
Green function is then written as 



Appendix B: Bulk Green's functions. 

If we set all the reflection amplitudes to zero in Eq. 
(|13p we obtain the bulk solution for an infinite graphene 
layer 

^bulk 



2hvc 



(Bl) 

Thus, the local propagator for a bulk of graphene is 



G'^""^(a; = x') 



2hv COS Q 

[cos^"' a(To + tanaCTy] 



G4, "-^ ^ ) " 2?w U^os <ao + tan<crj,} 



sS.bulk 
1 



{Etq + ATy + T, 



1 



+ (cos a,j(7o + tan af^dy) (g) — (i?fo + Aiy - 

(B3) 

In the heavily doped limit, the Green's function loses 
all structure in sublattice space, 



hvfi 



{Eto + Aiy) , (B4) 



(B2) and is equivalent to the bulk BCS Green function. 
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